rm(list=ls())
#library(readstata13)
library(foreign)
library(Zelig)
library(memisc)


dta <- read.dta("~/Dropbox/Chile Speeches/speeches_0607_c6.dta")

mod1 <- zelig(hi~ pdis+ lndistkm + keycomms + notrunagain + leadership + intralistmargin2 + marginlist2+freshmen+govt+died, model="negbin", data=dta) 

mod2 <- zelig(od~ pdis+ lndistkm + keycomms + notrunagain + leadership + intralistmargin2 + marginlist2+freshmen+govt+died, model="negbin", data=dta) 

mod1$setx(pdis=1.32,freshmen=0,died=0,govt=0,notrunagain=0,keycomms=0,leadership=0)
mod1$setx1(pdis=20,freshmen=0,died=0,govt=0,notrunagain=0,keycomms=0,leadership=0)

mod2$setx(pdis=1.32,freshmen=0,died=0,govt=0,notrunagain=0,keycomms=0,leadership=0)
mod2$setx1(pdis=20,freshmen=0,died=0,govt=0,notrunagain=0,keycomms=0,leadership=0)

mod1$sim(5000)
mod2$sim(5000)

m1fd <- mod1$getqi(qi="fd",xvalue="x1")
m2fd <- mod2$getqi(qi="fd",xvalue="x1")




par(mfrow=c(2,1),
	oma=c(2,2,0,2)+0.1,
	mar=c(2,0,2,0)+0.1,
	bty="n"
	)	
plot(density(m1fd),xlim=c(-50,150),ylim=c(0,0.025),yaxt='n',ann=F)
polygon(density(m1fd),col="gray",border="black")
abline(v=quantile(m1fd,c(0.05, 0.95)),lty=3)
abline(v=quantile(m1fd,c(0.5)),lty=5)
text(100,0.015,"Hora de Incidentes")

plot(density(m2fd),xlim=c(-50,150),ylim=c(0,0.035),yaxt='n',ann=F)
polygon(density(m2fd),col="gray",border="black")
abline(v=quantile(m2fd,c(0.05, 0.95)),lty=3)
abline(v=quantile(m2fd,c(0.5)),lty=5)
text(100,0.015,"Orden del Dia")

dta <- dta[dta$died!=1,]

mod.re <- zelig(retcongress ~ hi+ pdis+ lndistkm + keycomms + leadership + intralistmargin2 + marginlist2+freshmen+govt, model="logit", data = dta) 

mod.re$setx(hi = 5 , pdis = mean(dta$pdis,na.rm=T), lndistkm =  mean(dta$lndistkm,na.rm=T), intralistmargin2 =  mean(dta$intralistmargin2,na.rm=T), marginlist2= mean(dta$marginlist2,na.rm=T), freshmen=0,govt=0,keycomms=0,leadership=0)

mod.re$setx1(hi = 23 ,  pdis = mean(dta$pdis,na.rm=T), lndistkm =  mean(dta$lndistkm,na.rm=T), intralistmargin2 =  mean(dta$intralistmargin2,na.rm=T), marginlist2= mean(dta$marginlist2,na.rm=T),freshmen=0,govt=0,keycomms=0,leadership=0)

mod.re$sim(5000)

modre.fd <- mod.re$getqi(qi="fd",xvalue="x1")

modre.ev <- mod.re$getqi(qi="ev",xvalue="x1")

quantile(modre.fd,c(0.025,0.5,0.975))
quantile(modre.ev,c(0.025,0.5,0.975))



